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Numerical Computation of Fuel-Optimal, Low- and 
Medium- Thrust Orbit Transfers in 
Large Numbers of Burns 


Abstract 

This report presents two numerical methods considered for the computation of 
fuel-optimal, low-thrust orbit transfers in large numbers of bums. The origins of these 
methods are observations made with the extremal solutions of transfers in small numbers 
of bums; there seems to exist a trend such that the longer the time allowed to perform an 
optimal transfer the less fuel that is used. These longer transfers are obviously of interest 
since they require a motor of low thrust; however, we also find a trend that the longer the 
time allowed to perform the optimal transfer the more bums are required to satisfy 
optimality. Unfortunately, this usually increases the difficulty of computation. 

Both of the methods described use small-numbered bum solutions to determine 
solutions in large numbers of bums. One method is a homotopy method that corrects for 
problems that arise when a solution requires a new bum or coast arc for optimality. The 
other method is to simply patch together long transfers from smaller ones. An orbit 
correction problem is solved to develop this method. This method may also lead to a 
good guidance law for transfer orbits with long transfer times. 


I. INTRODUCTION 

Electric propulsion, with its high specific impulse, promises very low fuel 
consumption but it produces less thrust than its counterparts. If one wants to use electric 
propulsion, one needs to be prepared to tolerate the long transfer times that will be 
incurred. The greater time spent thru sting must be spent wisely if fuel savings are to 
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realized. Furthermore, the effects of Earth’s oblateness and atmospheric drag become 
more significant on the orbits of long transfer times. 

To spend the thrusting time wisely, we form an optimal control problem to 
maximize the mass at the end of the transfer. This, therefore, is our cost function 

J = m(t f ) ( 1 ) 

subject to the boundary conditions 

v( r (0), v(0),r(ry),v(/y)j = 0 (2) 


and the state dynamics 


f = v 



T 

m = 

8 o ^ sp 


(3) 

(4) 

(5) 


where e T is the thrust direction, a unit vector, and the thrust magnitude, T, is limited 
between zero and some maximum value T^, (i is the gravitational constant, g Q is the 
gravitational acceleration at sea-level, and / is the specific impulse of the motor. 
Sometimes gj sp is referred to as the exit velocity of the motor. If the boundary 
conditions referred to in Eqn. (2) are designed for the rendezvous problem, this results in 
the well-known bang-bang optimal control problem, discussed in detail by Lawden 1 . 
However, herein the boundary conditions are designed such that the initial and final 
points lie on the desired orbits without specifying the position, or true anomaly, on either 
orbit. 

As a brief review, the optimal thrust direction for this problem is 



where is found from the following differential equations 
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(7) 


K = m 


7" 


v -3 


(^ v T r)r 
— 3 — 


Xy — Xf 


a _ T » j _ T i» 

A m - — A v e T - ~ y| A > 


m 


m 


( 8 ) 

(9) 


The optimal thrust magnitude for this problem is a bang-bang solution. This is 
determined by applying the following switching law, Eqn. (10), to the switching function, 
Eqn.(ll). 


H s >0, T - 
H s < 0, T = 0 


( 10 ) 


// =fcvLja«_ 

™ 8o ! sp 


( 11 ) 


We are interested in solutions of this problem with long transfer times and, 
therefore, large numbers of bums. There are many methods that have been successively 
used to compute n-burn transfers, where n is anywhere from 1 to about 6. However, 
fewer methods successively compute transfers for larger values of n. These methods for 
the former attempt to solve the optimal control problem either directly, indirectly, or with 
a hybrid of the two. In this report, we will assume that a mostly indirect method, such as 
BOUNDSCO or MBCM or that of Brusch 2 , et. al, or of Redding 3 is being used. 

One idea to obtain interesting solutions is to first compute some n-bum transfer, 
where n is generally less than the number of bums initially desired. Using this as a 
starting point, increase the allowed transfer time and compute the new transfer. 
Obviously, it is expected that the desired transfer is relatively similar to the starting 
transfer. This homotopy method seems to work well as long as the number of bums 
performed in the transfer do not need to increase so that optimality is satisfied. For 
example, in many cases BOUNDSCO is unable to find a three bum solution when the 
two burn solution to an almost identical problem is given as the initial guess. The 
Direction Correction Method has been developed to attempt to alleviate this difficulty. 
It’s purpose is to find an n bum solution to an orbit transfer problem with allowed 
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transfer time ty + <5iy using an n-1 bum solution to the same problem but with allowed 
transfer time t^. 

Another idea is to patch together a set of n-bum transfers, where n is a small 
integer, usually unity, to produce an m-bum transfer, where m is the desired number of 
bums. This method requires that the sequence of transfer orbits be either guessed and 
iterated upon for optimality, or simply prespecified. From the theory of optimal control, 
this patched solution will be a suboptimal solution. However, possible analytical 
solutions for the one bum solution of two very close orbits may give a feedback guidance 
law. Since the drag model is only approximate for large numbers of bums it may be 
more important to have a good guidance law in terms of fuel-savings. 

H. Direction Correction Method 

The first idea, referred to herein as the Direction Correction Method, is based on 
the common homotopy strategy. A homotopy method, though slow in producing results, 
would be considered effective here as long as the number of bums does not change. It is 
expected, however, that one is going to be using this method to increase the transfer time 
so that the fuel consumed will be less. To understand the ensuing difficulty, we must 
study the history of a successful implementation of this homotopy method. 

All parameters describing transfers in this section and below have been 
nondimensionalized such that the gravitational constant, /t, is unity. This 
nondimensionalization is accomplished through two parameters, r* and m* with units 
of length and mass, respectively. These are chosen appropriately to the problem and may 
be, for example, initial semimajor axis and initial mass, respectively. The following 
equations detail the calculation of nondimensional parameters, denoted by the symbol, 
describing the transfer: 



(Usp)-sohpVr^ 


i - tf 


(12a) 

(12b) 

(12c) 
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The optimal transfer we will examine is a planar transfer under ideal gravity 
conditions. The transfer leaves an initial orbit with a semimajor axis of 2.239, 
eccentricity of 0. 1160, and an argument of perigee of -85.94°. The orbit to be entered has 
a semimajor axis of 7.000, eccentricity of 0.7332, and an argument of perigee of 1 14.6°. 
The motor used to perform the transfer delivers a thrust of 0.01386 with an exit velocity 
of 0.3898. The allowed transfer time is 73.33. This transfer performed in two bums is 
shown in Figure 1 with its corresponding parameters in Table I. It was computed using 
the multiple-shooting method of BOUNDSCO 4 . The switching function for this transfer 
is shown in Figure 2a. 



x 


Figure 1. Transfer in Two Burns for Burn Addition Demonstration 
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Table L Parameters of the transfer shown in Figure 1 
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Figure 2b Switching Function for a Two or Three Burn Transfer 
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The initial mass of the spacecraft was 1.6, the final mass is 0.5545. Now, suppose that a 
greater fuel savings is desired. As the allowed transfer time is increased from 73.33 to 
77.48 and then to 85.00, the shown sequence of switching functions (//y(f) in Figs. 2a-c) 
will result. These show a clear indication of a new bum/coast being anticipated in the 
optimal solution. The orbit transfer corresponding to the switching function in Fig. 2c is 
plotted in Figure 3. The parameters of this transfer are identical to that of Fig. 1 except 
that the transfer time is longer, tf= 85, see Table II for the listing. Also, note that the final 
mass of this longer transfer is indeed larger than the shorter transfer, indicating a greater 
fuel savings. 



x 


Figure 3 Transfer in Three Burns for Burn Addition Demonstration 
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Table n. Parameters of the transfer shown in Figure 3 
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It has been seen in many cases that local minima and maxima of the switching 
function will move down or up on the graph as we examine successive solutions. As in 
Fig. 2b, once this critical point becomes a root of the switching function, we reach a point 
where the number of bum/coasts is somewhat indeterminate. Is this, in Fig. 2b, a two- or 
three- bum extremal? There are only two bums of finite length but there is a third that is 
infinitely small. This indeterminacy shows itself as a discontinuity in the slope of a plot 
of the initial guess versus the homotopy variable, transfer time. Figure 4. 
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Figure 4 Plot of Initial and Final True Anomaly Values of Successive Solutions 
which Differ only in Transfer Time, tf. 


Figure 4 shows the initial and final true anomalies as a function of the allowed 
transfer time. The feature of interest here is the slope discontinuity (note that there is no 
point discontinuity) in both curves. The effect is not as prominent for the initial true 
anomaly as it is for the final, but it is still noticeable. As a result of this discontinuity 
there is difficulty in the homotopy method: the next solution may not converge because 
the method being used, based on the linear slope of previous points, is not calculating the 
correct initial state. To overcome this difficulty we must be able to compute the correct 
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slope, which should be the slope after j / = 77.48, so that the homotopy method can 
continue. 

The change in the initial state needs to be computed such that satisfaction of the 
boundary conditions is maintained and optimality is preserved. This problem shall be 
approached for the following general Two-Point-Boundary-Value-Problem (TPB VP): 

c( z (o)) = o m equations (13a) 

■#/)) = 0 m-n equations (13b) 

z(/) = f(z(/)) n equations (13c) 


where z(r) is the state consisting of the original state plus the Lagrange multipliers, f(/) is 

the right-hand side of the original state dynamics plus the Euler-Lagrange equations, and 
C(z(0)) and D(z(//)) are the boundary conditions for the initial and final orbits, 

respectively. 

Now, since we are interested in maintaining the boundary conditions, we set their 
variations equal to zero. First, the initial conditions from Eqn. (13a): 

= ^ <5z(0) = 0 (14) 

z(0) 


Next, a similar operation is performed on a vector describing the final conditions from 
Eqn. (13b). However, so that the initial state is referred to, it is necessary to invoke the 
transition matrix. 


d D = 


dD 

dz 

cT) 
dz . 

dz . 


in) 

<*/) 


dz[tf} = 0 

(dz(t f )+i(t f )dt f ) 
<D(0,r jr )fe(0)+^| i(t f )dt f = 0 

i { /) 


(15a) 

(15b) 

(15c) 


Here, d(-) denotes a variation with variable time and 0(0,//) is defined as the transition 
matrix, initialized at r=0, and evaluated at t=tf, where 
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(16a) 








Iz(r) 




(16b) 


Now at each switching time, fj (i=0,l...<7), the switching function must be satisfied. So, 
we set the variation of the switching function, H s ( z), equal to zero at each switching time. 

giving q scalar equations: 


II 

o 

II 

(17a) 

?§ * 
II 

(4>(0,i i )&(0)+i(f l > i ) = 0 

z(0 

(17b) 

II 

Dlj 

<f(0,/,)&(0) + ^ £ - i(t, >i = 0 

z(f:l *(*,) 

(17c) 


Consideration of the switching function also calls attention to a necessary correction in 
the transition matrix calculation. At each switching point, there is a discontinuity in 
f(z(t)) due to the thrust being turned on or off. This discontinuity results in a ‘jump’ term 
for 0(0, tf). To calculate this term, we again must set the total change in H s equal to zero. 


«,(«(%))■ 0 


dH s 



dz(r,) = 0 


(18a) 

(18b) 


Now, recognize that the total variation in the state at the switching time fj must be the 
same looking from either direction. This situation is illustrated in Figure 5. 


M*i) = 5z(t,- ) + z(f f )dt i 
= 5z(r i + ) + z(r, + )dr i 


(19a) 

(19b) 
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Figure 5. Illustration of Equations 19a and 19b 

Substitute Eqn. (19b) into Eqn. (18b) 

dH s =^- (&(,,-)+ z(< ->*,) = 0 (20) 


Equation (13) can be solved immediately for dr, which is then substituted into Eqns. (19a- 
b). This can manipulated to produce 

dH s 

dz 

dz 


5z(r, + ) = 8z(t~) + (z(r, + )- z(rf)) 



Equation (21) can be rewritten by inspection in terms of the transition matrix: 
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This is the jump matrix across the switching point 

We must recognize that these variations are considered in a range of transfer times 
across which the number of switching points changes. Specifically, this is an addition of 
a bum or coast arc. The situation is illustrated in Figure 6. The assumed change in the 
switching function is shown at the top of the figure. The nominal solution’s switching 
function has a touch point at tc=t a =t b . The solution with a slight different transfer time 
has two new switching points, t a + d t a and t b + dtb- The assumed change in one element 
of the state vector is shown at the bottom of Figure 6. The derivative, z(t), is assumed 
equal before and after the new addition and to the nominal value, z(t c ) . The slope during 

the new bum is denoted c. To relate the two solutions across the arc, we write the 
following equation. 

8z(t b + dt b )= 8z(t a + dt a ) + (c- z (t b ))(dt b - dt a ) (23) 

This relation has been verified using data from the example presented here. 
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Now, a model is required to locate the new switching points. We have looked 
into different models for this. The first model is a simple variational model, but unlike 
Eq (17a), second-order terms are considered. The equation of this model is 


A H s - -d 2 H s = -dz(t) T ^-^- 
s 2 * 2 Ka> dzdz 


Wa) 


dz{t a ) = 0 


= jH) + ^aK) T 


d 2 H< 


dzdz 


Hu) 


(Sz(t a ) + z(t a )dt a ) 




dzdz 


ko 


dz(t a ) + 25z(t a ) T ^-^ 


dzdz 


ko 


H‘a)dt a 


./ \jd 2 H s 
+z{t a ) 1 * 


dzdz 


HU) 




(24a) 

(24b) 


(24c) 


where the lesser of the two solutions is d t a leaving the other to be dr*. Unfortunately this 
model does not result in a sufficiently accurate answer for d t a and d /*. 

We have also attempted to model the situation through the information on the 
placement of t c + d t c . Since this point can be defined as the point of zero slope, we can 
find with an analog of Eq. (17). The solution is, therefore. 


dH t 


dt c = — 


dz 


k*e) 


MO 


dH< 


dz 


Hu) 


Hu) 


(25) 


To complete the model we need to have a point on the graph of A Hs and we need the 
curvature of Hs- The former can be had by rewriting Eq (17) for t c and evaluating it at 
d t c . We assume that the latter is well represented through a curve fit to the original 
switching function in the neighborhood of Hs(t c ), denoted by k. in the following equation 
for AH S . 


A H,~k(t-t c ) 2 + 




\ 

dz{t c ) 

) 


(26) 
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The solutions we are interested in are 


dt a = dt c - 


dt b = dt c + 


dH. 

s 

'N 

M'c) 
*(0 J 

k 

dH s 

dz 

\ 

M*c) 
<‘c) J 


k 


(27a) 


(27b) 


We have found that the solutions with this model are better than that with the previous 
model, but still not very accurate with errors greater than 10%. However, this accuracy 
may still prove to be well enough for BOUNDSCO to produce solutions. The intention 
here is merely to provide the TPBVP solver an initial guess closer to the n + 1 bum 
solution. 

Taking all of this together, a system of linear and non-linear equations can be 
written, starting with Eqs (14) and (15) 


dD 

dz 


dC 

dz 


\z(t f ) 


<&(0,t a +dt a )Sz(t a + dt a ) = 0 

z(0) 

<D(t* + dt b> t f )f(Sz(t a + dt a )) = - 




(28a) 

(28b) 


where f(<5z) refers to the right-hand side of Eq. (23) as a function of dz(t a + d t a ). The 
solution to this system is dz(t a + df a ), The transition matrix can be used to give the 
change in the initial state required to produce the desired solution. Then the variation of 
each switching point can be found one at time using Eqn. (17c). 

The solution information can easily be put into a form useful for a variety of 
numerical methods. For example, the change fiz(0) can be propagated through the 
transition matrix to calculate the changes at each node point for a multiple point shooting 
method. This method is still under development but shows promise as relatively simple 
way of getting to the n + 1 bum solution in the right direction. 

Once we have the ability to find optimal solutions with successively increasing 
transfer times, there is another characteristic of the extremals that may be encountered. 
Experience has shown that the length of the new bum will increase monotonically as the 
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transfer time is increased and usually the situation detailed above will be repeated so that 
the number of bums will increase again. However, there are cases where the cycle ends 
and the transversality condition, giving the optimal transfer time, is satisfied and there 
may be no nearby solution that has better performance. The following solution is an 
example. It is a descent trajectory from an orbit with a semimajor axis of 3.847, 
eccentricity of 0.02378, and an argument of perigee of 0°. The transfer terminates at an 
orbit with a semimajor axis of 1.500, eccentricity of 0.3333, and an argument of perigee 
of 0°. The motor used to perform the transfer delivers a thrust of 0.03 with an exit 
velocity of 1.313. The allowed transfer time is 19.05. This transfer performed in two 
bums is shown in Figure 7a. It also was computed using the multiple- shooting method of 
BOUNDS CO. The switching function for this solution is shown in Figure 7b. 



X 

Figure 7a: Two Burn Extremal with Transversality Converged 
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Table HI. Parameters of the Transfer Shown in Figure 7 
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Figure 7b: Switching Function for Two Burn Extremal in Figure 7a 

This solution was presented previously 5 , however, with one difference, oblateness and 
atmospheric drag were included in the dynamics. It was found that with these terms 
removed, the transversality condition could be converged. It was also observed that the 
initial and final points of the switching function were driven to zero. There is certainly 
no conflict here in terms of optimality: the initial and final points are now switching 
points. 

m. patched transfer method 


The second idea was inspired in part by the work of others. Zondervan, et. al 
made some simple guidance observations 6 , specifically that in some regions the primer 
vector is relatively constant in a velocity-fixed reference frame. This implies that a 
simple control law is available in some cases. Marec presents a solution to the orbit 
correction problem 7 . This motivated a notion that solutions to linearized and/or 
approximated problems were available. In this spirit a solution was obtained for the 
optimal transfer between two close orbits. The transfer leaves a circular initial orbit with 
a radius of 1.038. The orbit to be entered has a semimajor axis of 1.069, eccentricity of 
0.02633, and an argument of perigee of -50°. The motor used to perform the transfer 
delivers a thrust of 0.01438 with an exit velocity of 0.3861. The allowed transfer time is 
1.553. This transfer is performed in one bum and is shown in Figure 8a. It was 
computed using the multiple- shooting method of BOUNDSCO. The switching function 
for this transfer is shown in Figure 8b. 
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Figure 8a: One Burn Transfer Between Close Orbits: An Example of a Solution 

with a Simple Optimal Control 
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Table IV. Parameters of the Transfer Shown in Figure 8a 

Most interesting about this transfer is the simplicity of the control. Over this short 
transfer between a circular orbit and a close target orbit, the optimal control of the thrust 
angle is linear in time. And, in addition, we find that the control direction is almost 
coincident with the velocity direction. 
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Figure 8b: Switching Function for One Burn Transfer Shown in Figure 8a 



Figure 8: Plot of Thrust Direction, the Optimal Control, Alongside the Angle of 

the Velocity Vector. 

To match this transfer analytically, a modified optimal control problem is 
considered. The dynamics for this problem are again the equations of orbital motion, 
however, this time the state is defined relative to the initial orbit. Assuming that the 
distance from the initial orbit is small compared to the radius of the initial orbit, we 
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ignore all terms to the order of (<5r/p) 2 . This assumption results in the following 
dynamics: 


$ 

II 

* 

(29a) 

T u(5r*p) p _ 

<5v = — e r + 3 — — f— ^ p--^<5r 
m p 3 p 

(29b) 

T 

m = 

So^sp 

(29c) 

Here, 5r=[jc y] T and 5v=[u v] T , e T is the thrust direction, T is the thrust, m is the mass, p 
is the gravitational constant, and p represents the initial orbit which satisfies identical 
dynamics but without the thrust term. Now, assuming that the initial orbit is circular, 

these can be rewritten as: 


x = u 

(30a) 

y = v 

(30b) 

u = —e x + -^-[3 (xcos(at) + ysin((ot))cos((dt) - x] 

(30c) 

v = —e y + [3(;t cos(ox) + y sin(a)t)) cos(eot) - y] 

(30d) 

T 

m = 

So^sp 

(30e) 


Writing the Hamiltonian for this system and evaluating the Euler-Lagrange equations 
results in the following differential equations involving the costates: 


H 

1 

II 

(31a) 

1 

II 

> 

(31b) 

X x = - 3(3 cos 2 (cot) - 1) - A v 3 cos{(Ot)sin{(Ot) 

(31c) 

X y = -X u -^j(3cos((ot)sin(ax))- X v -^j{3sin 2 (cot)-lj 

Old) 

i - T 

m 2 ^X 2 + X v 2 

Ole) 


We also learn that the control, e T is 
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(32) 


e T - 


1 

PM 

-\j X 2 + X v 2 

kJ 


and the control T is bang-bang, governed by the switching function, Hj, as 
_ V^« 2+ ^v 2 X m 

Ll'Y — 

® 8cJsp 

Hj> 0, T = 

h t < o, r = o 


(33) 

(34) 


Pleasantly, Eqns. (31) happen to be the equations for the costates on a coast arc 
coinciding with the initial orbit. In fact, this result is not limited to the assumption of a 
circular orbit. The coast arc costates have been solved by Lawden and other authors 8 * 9 . 
It also can be shown that Eqns. (31) are, in fact, identical to Eqns. (30), without the thrust 
terms, up to sign. Therefore, once we solve the system in Eqns. (31) we have the 
homogeneous solution to the system in Equations (30). Now to solve the differential 
Eqns (26), they must first be rewritten in a more useful form: 




= - 3 1 


cos(cot ) 
sin(cot) 


[cos(cot) 



(35) 


where Now, define vectors e p (f) and e^r), as the radial and circumferential 

directions associated with the initial orbit over time t. This can now be written as 

X = 3le p e p r X-!X (36) 

T 

where X = [X u X v ] . Multiply Eqn. (36) by e p T and e<y T , respectively to obtain 

e p T X = 3 le p T X - le p T X = 2 le p T X = 2 co 2 e p T X 
eJX = 3leJe p e p T X-ltJX = -le (0 T X = -CD 2 e 0) T X 


(37a) 

(37b) 


To complete the simplifications, it is necessary to obtain an expression for left-hand side 
of Eqn. (36) in terms of e p and . That expression is 
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(38) 


= |i!p — 2<t)A ei j - ft) Xp jcp + ^A^, + 2ft)Ap — ft) A(0 jcju 

Using this expression, Eqns. (37) become 

Ap -2 ci)i a - ci) 2 X p = 2co 2 X p 
A© + 2 (oi p - ft ) 2 A c , = -fl> 2 A» 


This can be represented with a matrix differential equation. 


V 


■ 0 

0 

1 

0 ' 

V 



0 

0 

0 

1 

^Q) 

K 


3ft) 2 

0 

0 

2 ft) 

Ai 

A 2 _ 


0 

0 

-2 ft) 

0 

A2 


where Aj= dkp/dt and X^= dkjdt. The solution to this system is 


'V 


' 2 ' 


cos(cot) 


sin(ox ) 


O' 

^0) 

-a 

-Scot 

+ Z> 

-2sin( cot) 

+ c 

2 cos{(Ot ) 

+ d 

1 

Ai 

0 

-Q)sin(i ot) 

cocos{(Ot) 

0 

_a 2 _ 


3ft) 


-2cocos(cot) 


-2cosin(a)t) 


0 


(39a) 

(39b) 


(40) 


(41) 


where a, b, c, and d are independent constants. The vector A can be interpreted as the 
thrust direction in a reference frame fixed to the radius of the initial orbit, referred to here 
as the initial orbit reference frame. From the solution above, Eqn. (41), we see that there 
are four modes of the thrust direction. The mode associated with d is fixed with respect 
to the initial orbit reference frame. The mode associated with a is not fixed to that frame 
but is very simply described in it. The last two modes do not seem as well described in 
this frame. 

To be sure, we would like to see that the approximate state dynamics given in 
Eqns (29) and (30) closely match those given in Eqns (3,4,5). To validate the 
approximate dynamics, it was simplest to simulate both systems using the same control. 
The most obvious choice for this control is the optimal control from the transfer in Fig. 
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7a. Figure 9 shows the results of the simulation. In this figure, “Delta-” states refer to 
the states from Fig. 7a with the initial orbit states subtracted, producing the desired plot 
for &*. The “X1,Y1,” etc. states refer to the states obtained by integrating Eqns. (30). 
The results seen in this figure are very promising: there is almost exact agreement 
between the two state histories. In fact, the worst error between the two at the end of the 
transfer is only about 1.5%. 


■s — Delta-X — 6— Delta-U — XI — *— U1 

e— Delta-Y — v— Delta-V — *— Y1 — VI 



Figure 9: Validation Plot for the Dynamics in Equations (29) for the Transfer 

shown in Figure 7a 

IV. Conclusions 

The development of the Direction Correction Method is proceeding rather well. 
The ideas that it is based upon have been validated individually. At this point, the only 
weak link is the prediction of the new switching points. Testing of the method will be 
required in order to determine just how critical is the accuracy of that prediction. 
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The Patched Transfer Method is very promising. The dynamics resulting from 
assumptions made closely matches the dynamics before the approximations. Also, the 
simplicity of the resulting optimal control problem promises a state feedback guidance 
law. The usefulness of these results will outweigh the loss in accuracy. However, much 
more analysis must be performed to completely validate the linearized problem and its 
solution. Specifically, the approximate optimal control solution must be compared to 
exact solution; based on the agreement of the state, positive results are expected, but they 
must be verified. 
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